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1. Introduction 



The atomic nuclei have been historically treated as collections of protons and neutrons. The 
great success of the nuclear shell model since 1949 [|l|, |2|], explaining the nuclear magic numbers 
and detailed spectroscopy, has established that protons and neutrons are very good effective degrees 
of freedom at the nuclear energy scale of a few MeV. Nonetheless, 60 years later, we know for 
certain that protons and neutrons are made of quarks and gluons whose laws are governed by QCD. 
It is a great challenge to quantitatively understand the structure and property of known nuclei based 
on the first principle of QCD. This direct approach will be more important and indispensable if we 
are to extract reliable predictions for experimentally unknown nuclei in the neutron rich regions 
of the nuclear chart. In this article we address the fundamental question in the research in this 
direction, namely the binding energies of nuclei. 

Lattice QCD study of multi-baryon states goes back a long time, starting with H dibaryon 
in the 80's |Q] and deuteron in the early 90's [Q]. More recently, exploration of three baryon states 
began [|6|]. So far, however, there is no established evidence supporting bound state formation in 
these channels. An exception is a model study in the strong coupling limit of lattice QCD [^]. 

We attempt to go a step further in mass number and examine the helium nuclei, especially 
^He with the mass number A = 4. Besides the obvious physical interest as the first natural element 
beyond hydrogen, it is also the system where technical difficulties of fermion contractions specific 
to nuclei with a large mass number appear in a non-trivial way. On the other hand, the binding 
energy drops down to a large value of AE = 28.3 MeV for the ^He nucleus, making us hopeful that 
observing the bound state nature might be easier than the lighter nuclei. 

The organization of this article is as follows. In Sec. § we review previous studies for bound 
states in multi-baryon systems from lattice QCD. The computational issues with studies of multi- 
baryon states and their solutions employed in this work are briefly explained in Sec. ^ The simula- 
tion details and the results for the ^He and ^He channels are presented in Sec. Q A brief summary 
and a look toward future are given in Sec. |[ The results in this article have been reported in 
Ref. [§. 



2. Historical perspective 

Bound states in multi-baryon systems have been investigated by several studies in lattice QCD. 
For systems with two baryons, the first study was the search for the H dibaryon. According to 
Jaffe the H dibaryon having strangeness S = —2 and isospin 7 = channel was expected to 
have a large binding energy of 0(100) MeV. Most of the quenched lattice QCD studies [0, ^, [l^. 



11, 12] concluded that the H dibaryon bound state does not exist. Recently NPLQCD Collaboration 



observed a small, negative energy shift, E/^ — 2m/^ = — 4.1(1. 2)(1.4) MeV [13], in this channel. 
They concluded, however, that the evidence is not strong enough to establish the existence of the 
H dibaryon, and that further study is necessary with different volumes. The latter point is related 
to the computational problem of the nucleus calculation, which we will discuss in the next section. 

Deuteron is a bound state of two nucleons in the and 7 = channel. Nucleon-nucleon scat- 
tering in this channel and also in the '^o channel was first studied in quenched QCD [S, 14]. This 



work was followed by a partially-quenched mixed action [15] and Nf = 2 + 1 anisotropic Wilson 
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Figure 1: Scattering length for ^Sq (left) and (right) channels. Circle, square, diamond, and triangle 
denote results for quenched mixed action jl^, two-nucleon wave function jl^, and anisotropic 

Wilson [pj| calculations, respectively. 



action 0130 simulations. Extraction of nuclear force between two nucleons has been investigated in 
quenched and 2+1 flavor QCD [16, [ij, 18]. Results for the scattering lengths oq from these studies 
are summarized in Fig. |[ The scattering lengths in the two channels are almost identical in each 
group. The results, however, have large discrepancies between the groups. An even more prob- 
lematic issue is that the absolute value of the lattice results is much smaller than the experimental 
values, qq = 23.7 fm and ao = —5.47 fm for the '^o and channels, respectively. The lattice 
results do not show strong dependence on the pion mass at the region where the calculations were 
carried out, ni;i:^0.3 GeV. In order to explain the experimental values, the scattering lengths have 
to vary significantly when calculations near the physical quark mass are carried out in future. We 
should also note that all these studies assumed that the deuteron state is not bound for the heavy 
pion mass employed in the calculations. 

Recently not only two-baryon systems but also three-baryon systems have been investigated 
using lattice QCD. NPLQCD Collaboration has tried a feasibihty study of three-baryon systems 
focusing on the 'ZP'EPn and the nnp (triton) channels. They found both interactions to be repul- 
sive ||6, 19], which indicates that the triton is not bound for the parameters taken for the calculation. 

In this conference several studies of two- and three-baryon systems were reported. HALQCD 
Collaboration studied (i) the energy dependence of the nuclear force [pO|], (ii) the nuclear force 
in the flavor SU(3) limit [^], (iii) extraction of the two-baryon forces in a coupled channel with 



the variational method [22, 23], and (iv) an exploratory study of extraction of the three-nucleon 



force [24]. In multi-meson systems, NPLQCD Collaboration proposed a recursion relation ap- 



proach for multi-meson correlation functions [ |25| , [26[ ] to largely reduce the computational cost of 
the correlation functions. 



3. Computational issues with nuclei 

There are several computational difficulties in the calculation of the multi-baryon bound state 
in lattice QCD. They are : 1) exponential increase of statistical error, 2) factorial growth of fermion 
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Wick contractions, and 3) identification of bound state. While we avoid the first one by an unphys- 
ical heavy quark mass, we propose solutions for the second and third problems. Let us discuss each 
in turn. 



3.1 Exponential increase of statistical error 

An estimate of the statistical noise to signal ratio for the correlation function of the nucleus 



consisting of Nn nucleons is known [27] to be proportional to 

1 r 3 " 



:exp Nn 



(3.1) 



where niji and mj^ are the masses of the pion and nucleon, respectively, A^meas is the number of 
measurement, and t is the separation between the source and sink time slices. The statistical error 
exponentially increases as the number of nucleon increases as well as when the quark mass de- 
creases. We aim to treat helium nuclei in this work, so that is fixed to four and three for ^He 
and ^He channels, respectively. Since our main aim is to explore nucleus calculations, and since 
the difficulty of controlling statistical fluctuations toward the region of lighter pion mass is well 
known, we use the heavy quark mass corresponding to nin = 0.8 GeV. Even then we had to carry 
out 0(10^^) measurements. 

While this strategy would be acceptable for a feasibility test of calculation of nucleus, we need 
novel methods to solve this problem for a more realistic calculation near the physical quark mass. 
We leave this task in future. 



3.2 Factorial growth of Wick contractions 

Another computational problem for multi-nucleon systems is a factorially large number of 
Wick contractions of quark-antiquark fields required for evaluations of the nucleus correlation func- 
tions. A naive counting would give {2Np +N„) ! {2N„ +Np) ! for a nucleus composed of Np protons 
and A'^ neutrons, which quickly becomes prohibitively large beyond three-nucleon systems, e.g., 
2880 for ^^He and 518400 for "^He. 

This number, however, contains equivalent contractions under the permutation symmetry in 
terms of the protons or the neutrons in the interpolating operator. We can reduce the computational 
cost by avoiding the redundancy. In the case of the ^He nucleus which consists of the same number 
of protons and neutrons, the isospin symmetry also helps us reduce the necessary contractions. 
After a scrutiny of the remaining equivalent contractions by a computer we find that only 1107 (93) 
contractions are required for the ^He (^He) nucleus correlation function. We have made a numerical 
test that the result with the reduced contractions reproduces the one with the full contractions on a 
configuration. 

For an additional technique to save the computational cost of the nucleus correlation functions, 
we make a block of three quark propagators where a nucleon operator with zero spatial momentum 
is constructed in the sink time slice. In this procedure we can incorporate the permutation symmetry 
of two up (down) quarks in a proton (neutron) sink operator. This is a simple trick to calculate 2^-^ 
contractions simultaneously. We also prepare several combinations of the two blocks which are 
useful for the construction of the nucleus correlators. 
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L 


A'^conf 


^meas 


accept. (%) 


[GeV] 


MN [GeV] 


24 


2500 


2 


93 


0.8000(3) 


1.619(2) 


48 


400 


12 


93 


0.7999(4) 


1.617(2) 


96 


200 


12 


68 


0.8002(3) 


1.617(2) 



Table 1: Number of configurations (A^conf)^ number of measurements on each configuration (Nmeas), accep- 
tance rate in the HMC algorithm, pion mass (niji) and nucleon mass (niN). 



3.3 Identification of bound state 

A general issue with numerical calculations for exploring bound state formation is to dis- 
tinguish the physical binding energy from the energy shift due to the finite volume effect in the 



attractive scattering system g9|, The problem is made more difficult for nuclei because 
the binding energy AE of the nucleus consisting of Nn nucleons with the mass is very tiny 
compared with the mass M of the nucleus: AE/M ~ O(10~^) with AE" = NNmN — M. 

One way to solve the problem is to investigate the volume dependence of the measured energy 
shift: In the attractive scattering system the energy shift is proportional to 1 /L^ at the leading order 



in the 1/L expansion [ pSj pl| ], while the physical binding energy remains at a finite value at the 
infinite spatial volume limit. 

If the volume is not large enough, it is difficult to distinguish a constant from a 1 /L? behavior 
in the energy shift. Thus, in our simulation we employ large volumes as much as possible, and 
choose three spatial extents corresponding to 3.1, 6.1 and 12.3 fm. The largest two volumes are 
much larger than those employed in current numerical simulations. They should provide sufficient 
room for the helium nuclei. 



4. A quenched calculation of Helium nuclei 

4.1 Simulation details 

We carry out calculations on quenched configurations generated with the Iwasaki gauge ac- 
tion [32] at j8 = 2.416 whose lattice spacing is a = 0.128 fm determined with ro = 0.49 fm as an 
input [Q. We employ the HMC algorithm with the Omelyan-Mryglod-FoUc integrator [|^, ^]. 
The step size is chosen to yield reasonable acceptance rate presented in Table [l[ We take three lat- 
tice sizes, xT = 24^ x 64, 48^ x 48 and 96^ x 48, to investigate the spatial volume dependence 
of the energy difference between the ground state of the nucleus channel and the free multi-nucleon 
states. The physical spatial extents are 3.1, 6.1 and 12.3 fm, respectively. 



We use the tadpole improved Wilson action with csw = 1-378 []33|]. As discussed in the pre- 
vious section, since it becomes harder to obtain a reasonable signal-to-noise ratio at lighter quark 
masses for the multi-nucleon system, we employ a heavy quark mass at K" = 0.13482 which gives 
ifijc = 0.8 GeV for the pion mass and =1.6 GeV for the nucleon mass. Statistics are increased 
by repeating the measurement of the nucleus correlation functions with the source points in differ- 
ent time slices on each configuration. The numbers for the configurations and the measurements 
on each configuration are summarized in Table We separate 100 trajectories between each mea- 
surement with T = 1 for the trajectory length. The errors are estimated by the jackknife analysis 
choosing 200 trajectories for the bin size. 
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The quark propagators are solved with the periodic boundary condition in all the spatial and 
temporal directions, and using the exponentially smeared source 

cj[{x,t)=Y^Ae-^\^-y\q{y,t) (4.1) 

y 

after the Coulomb gauge fixing, q is the quark field at the source time slice and A, B are the smearing 
parameters. On each volume we employ two sets of the smearing parameters: (A,B) = (0.5,0.5) 
and (0.5,0.1) for L = 24 and (0.5,0.5) and (1.0,0.4) for L = 48 and 96. Effective mass plots with 
different sources, which are shown later, help us confirm the ground state in the nucleus channel. 
Hereafter the first and the second smearing parameter sets are referred to as "Si .2", respectively. 

The interpolating operator for the proton is defined as pa = £abc{[ii-a\Cy5db)u'^ where C = 7472 
and a and a,b,c aie the Dirac index and the color indices, respectively. The neutron operator is 
obtained by replacing u" by in the proton operator. To save the computational cost we use the 
nonrelativistic quark operator, in which the Dirac index is restricted to upper two components. 

The "*He nucleus has zero total angular momentum and positive parity 7^ = 0+ with the isospin 
singlet 7 = 0. We employ the simplest '^He interpolating operator with the zero orbital angular 
momentum L = 0, and hence J = S with 5 being the total spin. Such an operator was already given 



long time ago in Ref. [p6[], 

'Ue = {xri-xri)/V2, (4.2) 

where 

X = ([+-+-] + [- + -+]-[+--+]-[- + +-])/2, (4.3) 
X = {[+-+-] + [- + -+] + [+ +] + [- + +-] - 2[+ + — ] - 2[- - ++])/^(4.4) 

with +/— being up/down spin of each nucleon. 17,7] are obtained by replacing +/— in by 
p/niox the isospin. Each nucleon in the sink operator is projected to have zero spatial momentum. 
We also calculate the correlation function of the ^He nucleus whose quantum numbers are 



= i^, 7 = i and 7^ = 4. We employ the interpolating operator in Ref. [37 1, 



2 ' ' ~ 2 — 2 • 

^He = {\p^n^p+) - \p+n+p^) + \n+p+p^) - \n+p^p+) + \p+p^n+) - \p^p+n+)) /Vg, (4.5) 
with the zero momentum projection on each nucleon in the sink operator. 
4.2 "^He channel 

Let us first present the results for the "^He channel. The left panel of figure § shows the effective 
mass plots of the "^He nucleus correlators with the Si 2 sources on the (6.1 fm)^ spatial volume. We 
find clear signals up to f « 12, beyond which statistical fluctuation dominates. The effective masses 
with the different sources show a reasonable agreement in the plateau region. The consistency is 
also shown in the exponential fit results in the plateau region as presented by the solid lines in the 
figure. 

In order to determine the energy shift AEl precisely, we define the ratio of the ^He nucleus 
correlation function divided by the fourth power of the nucleon correlation function. 
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Figure 2: Effective masses for ^He (left) and ^He (right) correlation functions with (circle) and ^2 (square) 
sources at spatial extent of 6.1 fm. Fit results with one standard deviation error band are expressed by solid 
lines. 



0.04 
0.02 


-0.02 
-0.04 
-0.06 



T 1 1 T 1 1 1 1 1 1 T— 



He 



< 1 



T 1 1 1 T— 



12 



16 



20 




Figure 3: Effective energy shifts for ^He (left) and ^He (right) channels in a convention of —A.E'f^ with Si 
(circle) and ^2 (square) sources at spatial extent of 6.1 fm. Square symbols are slightly shifted to positive 
direction in horizontal axis for clarity. Fit results with one standard deviation error band are expressed by 
solid lines. 



where G4He(0 ^"^^ G^it) are obtained with the same source. The effective energy shift is extracted 
as 

_A£f = lnf^^^y (4.7) 

once the ground states dominate in both of the correlators. In the left panel of Fig. || we present 
time dependence of for the 5i.2 sources, both of which show negative values beyond the 

error bars in the plateau region of 8 < f < 11. Note that this plateau region is reasonably consistent 
with that for the effective mass of the ^He nucleus correlators in the left panel of Fig. The signals 
of -d£f are lost beyond t ^ 12 because of the large fluctuations in the "^He nucleus correlator. 
We determine AEl by exponential fits of the ratios in the plateau region, f = 8 — 12 for Si and 
f = 7 — 12 for S2, respectively. We estimate a systematic error of AEl from the difference of the 
central values of the fit results with the minimum or maximum time slice changed by ±1. 

Table || summarizes the numerical values of the energy shift AEl at three spatial volumes, 
where the statistical and systematic errors are presented in the first and second parentheses, respec- 



7 



Calculation of Helium nuclei in quenched lattice QCD 




Figure 4: Spatial volume dependences for —AEi = M — NNniN in GeV units for ''^He (left) and ^He (right) 
nuclei with (open square) and S2 (open diamond) sources. Statistical and systematic errors are added 
in quadrature. Diamond symbols are slightly shifted to positive direction in horizontal axis for clarity. 
Extrapolated results to the infinite spatial volume limit (filled circle) and experimental values (star) are also 
presented. 



lively. The volume dependence of AEl is plotted as a function of in the left panel of Fig. ^. 
The errors in the figure are evaluated from the statistical and systematic errors added in quadrature. 
In the following discussions in this subsection we use the combined error. The results for the 5i,2 
sources are consistent within the error bars. We observe little volume dependence for AEl indicat- 
ing a bound state, rather than the dependence expected for a scattering state, for the ground 
state in the "^He channel. 

The physical binding energy AE defined in the infinite spatial volume limit is extracted by a 
simultaneous fit of the data for the S\2 sources employing a fit function of AE + C/L^ with AE 
and C free parameters. The term is added to allow for contamination of scattering states. A 
systematic error is estimated from the difference of the central values of the fit results using the data 
with the different fit ranges in the determination of AEl. The result for AE is 0.0180(62) in lattice 
units, which is 2.9 a away from zero as shown in the left panel of Fig. ^. We also try a pure bound 
state fit allowing for an exponentially small finite size correction: AE and AE + Cie^''^^ with AE 
and Ci 2 free parameters. We find all the results are in agreement with reasonable values of x^- 

Based on these analyses we conclude that the ground state of the measured four-nucleon sys- 
tem is bounded. An encouraging finding is that AE = 27.7(9.6) MeV with a^' = 1.54 GeV agrees 
with the experimental value of 28.3 MeV. However, we do not intend to stress the consistency be- 
cause our calculation is performed at the unphysically heavy pion mass, rrijc = 0.8 GeV, and the 
electromagnetic interactions and the isospin symmetry breaking effects are neglected. 

4.3 ^He channel 

The results of effective mass and the effective energy shift for the ^He channel with the Si 2 
sources are shown in the right panel of Figs. ^ and ^, respectively. The statistical error is slightly 
smaller than those for the '^He channel. We determine AEl for the ^He channel as in the "^He 
channel, whose results are presented in the right panel of Fig. ^ and Table ^ The trend of the 
volume dependence is similar to the '^He channel case. A simultaneous fit of the data for the 51.2 
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L 






[MeV] 






^He(5i) 


'He(52) 


'He(5i) 


^He(52) 


24 


28(14)(11) 


46.8(7.3)(1.6) 


19.0(6.3)(6.0) 


23.2(3.2)(0.5) 


48 


27(14)(05) 


36(12)(04) 


16.6(6.9)(3.2) 


19.5(5.6)(2.3) 


96 


24(18)(12) 


24(14)(03) 


19.0(7.6)(4.9) 


18.4(6.1X1.9) 


OO 


27.7(7.8)(5.5) 


18.2(3.5X2.9) 



Table 2: Energy shifts for ''^He and ^He channels on each spatial volume. Extrapolated results to the 
infinite spatial volume limit are also presented. The first and second errors are statistical and systematic, 
respectively. 



sources with a fit function of AE + C/Lr' yields a finite value of AE = 18.2(4.5) MeV, with the 
combined errors as in the "^He channel, in the infinite volume limit. This means the existence 
of a bound state in the ^He channel. Our result for AE, however, is about twice larger than the 
experimental value of 7.72 MeV. The main reason could be the heavy pion mass employed in this 
calculation. 

As an alternative way to view this result, we compare the binding energies normalized by the 
atomic number: AE/Nn = 6.9(2.4) MeV and 6.1(1.5) MeV for the "^He and ^He nuclei, respec- 
tively. At our unphysically heavy pion mass, the three and four nucleon systems do not show the 
experimental feature that the binding is stronger for "^He than for ^He. We consider that this is 
mainly caused by the heavy quark mass used in this calculation. 



5. Toward further progress 

We have addressed the issue of nuclear binding for the case of '^He and ^He nuclei. We have 
shown that the current computational techniques and resources allow us to tackle this issue. Albeit 
in quenched QCD and for unphysically heavy pion mass, we are able to extract evidence for the 
bound state nature of the ground state and the binding energies for these nuclei. 

It is encouraging that our results for the binding energies and the experimental values are 
of same order of magnitude. There are several issues which need clarification, however. Our 
conclusion for ^He seems to disagree with the recent result of NPLQCD Collaboration [^, p^ ]. 
While the two simulations have been carried out under different parameters, e.g., quark mass, 
number of flavors, and volumes, further work is needed to obtain a consensus in this chamiel. 
Furthermore we should gain a deeper understanding on how the helium nuclei forms a bound state 
at such a heavy quark mass. Study of the nuclear force extracted from the wave function and 
looking at the simplest nucleus, deuteron, might help to understand our results. 

A future direction of primary importance is to investigate the quark mass dependence of the 
binding energies of the nuclei. There are several model studies of the quark mass dependence of 



the nuclear binding energies [38] which suggest that the quark masses play an essential role in a 
quantitative understanding of the binding energies. 

Another important issue is development of a strategy to calculate nuclei with larger atomic 
numbers. The number of Wick contractions quickly diverges as the atomic number increases, 
even if the redundancies are removed with various symmetries and techniques employed in this 



9 



Calculation of Helium nuclei in quenched lattice QCD 



work. At this conference, NPLQCD Collaboration presented a set of recursion relations |25, 26] 
for correlation functions for the ?i-meson system, and also for more complex systems with multi- 
species, such as n-n and m-K systems. Similar recursion relations in multi-baryon systems might 
solve the problem. We leave it to future work. 
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